clear 
est clear

********************************************************************************
*	Programs
******************************************************************************** 

*** program for the OLS regressions... 
capture program drop regs
program define regs, eclass

	*** Syntax 
	syntax varlist(numeric) [, treatment(varlist numeric) rounds(numlist) ctrls(varlist numeric)]
 
	*** loop over outcomes and rounds
	qui foreach y of local varlist {
		foreach r in `rounds' {
		    
			*** create the block variables, always centred for sample with nonmissing outcome 
			capture drop blocks*
			qui tab block, gen(blocks)
			qui foreach v of varlist blocks* {
				replace `v' = . 			if missing(`y') & round ==`r'
				su `v' 						if round == `r' & !missing(`y')
				replace `v' = `v' - r(mean) if round == `r' & !missing(`y')
			}
			
			*** Model, includes: 
			*	- treatment,
			*	- block FE fully interacted,
			*	- clustered SE's at club level,
			*	- tracking weights (=1 in RR and 12m fups)
			* 	- controls fully interacted (when adjusted)
			local model_name = "`y'_`r'"
			* if (length("`ctrls'")>0) local model_name = "`y'_`r'_adj"
			
			eststo `model_name': reg `y' b0.`treatment'##c.(`ctrls' blocks*) [pw=w] if round==`r', vce(cluster club)
			
			*** save stats 
			* control mean and SD 
			qui su `y' [aw=w] if `treatment'==0 & round==`r' 
			estadd scalar c_mean = r(mean) 	
			estadd scalar c_sd = r(sd) 
			* test coeffs p-values, when treatments are combined
			if ("`treatment'"=="treat_ther") { 
				* p-value for IPT-G
				test 1.treat_ther == 0 
				estadd scalar pval_ther_any = r(p)
			}
			* p-value for the two treatments and comparison stats 
			else if ("`treatment'"=="treat") {
			    * p-value for IPT-G
				test 1.treat == 0 
				estadd scalar pval_ther = r(p)
				* p-value for IPT-G+
				test 2.treat == 0 
				estadd scalar pval_cash = r(p) 
				* diff between treaments 
				lincom 2.treat - 1.treat
				estadd scalar est_diff = r(estimate) 
				estadd scalar pval_diff = r(p) 
			}
		}
	}
	drop blocks* 
end


*** program to extract q values for the outcomes and treatment specified, within each round
capture program drop get_qvals
program define get_qvals, eclass

	*** Syntax 
	syntax varlist(numeric) [, rounds(numlist) tname_in(string) tsuff_out(string) method(string)]
	
	*** Create vector of p-values
	* matrix to store them
	local varcount : word count `varlist' 
	matrix pvals_tab = J(`varcount',1,.)  
	matrix colnames pvals_tab = pvals
	
	* loop over outcomes and rounds
	qui foreach r in `rounds' {
	    local i = 0
		foreach y of local varlist {
		    local ++i
			* restore estimates
			est restore `y'_`r'
			est replay
			* collect pvals in matrix 
			if ("`tname_in'" == "est_diff") {
				lincom 2.treat - 1.treat 
				matrix pvals_tab[`i',1] = `r(p)'
				matrix list pvals_tab				
			}
			else {
				matrix list r(table)
				matrix pvals_tab[`i',1] = r(table)["pvalue","`tname_in'"]
				matrix list pvals_tab			    
			}
		}
		* save in a dataset 
		capture drop pvals`r'
		svmat pvals_tab, names(col)
		rename pvals pvals`r'
		
		*** get q-values  
		capture drop qvals`r'
		qqvalue pvals`r', method(`method') qvalue(qvals`r') 
	}
	
	*** add to estimates 
	* loop over outcomes and rounds
	qui foreach r in `rounds' {
	    local i = 0
		foreach y of local varlist {
		    local ++i
			* restore estimates
			est restore `y'_`r' 
			* add as a scalar
			local qv = qvals`r' in `i'
			di `qv'
			estadd scalar qval_`tsuff_out' = `qv'
		}
	} 
end






********************************************************************************
*	I - Impacts of Therapy vs. Control
********************************************************************************

use "${data}/Analysis_wide.dta", clear
distinct block club

*** is it easier in long mode? I think it might be...
keep         age_current? w? *_b block cr_id club treat treat_ther treat_cash
reshape long age_current  w, i(block cr_id club treat) j(round)

*** Merge in relevant covid vars 
merge 1:1 cr_id round using  "$data\SMU_CleanData_Panel_Public", nogen assert(match) keepusing(ct_* c19_*) // keep(match) 


**** Some Descriptives of the CT ***********************************************

*** Do they still have any savings?
* only one third... 
tab ct_anyleft round, col 

*** How many gave cash to family/friends? 
* most had to give some
tab ct_hadtogive round , col

*** Could they choose how to spend it? 
* most still got to choose how to spent 'some' of it...
tab ct_choosespent round , col

*** Thinks Covid affected how to spend 
tab ct_covideffect round, col

*** Thinks the CT helped 
tab ct_helpcovid round , col
* and how? 
foreach v of varlist ct_helpC19_buymeds ct_helpC19_buyfood ct_helpC19_payrentbills ct_helpC19_supportfam ct_helpC19_invest ct_helpC19_other {
	tab `v' round , col
}


**** Impact of Covid ***********************************************************

*** prep covid vars 
* rename c19_* *
foreach v of varlist c19_getangry c19_argue c19_talkthrough c19_helpchores c19_fear {
	replace `v'=. if `v'<0
	tab `v' round , col 
}
* increased stress agree/partially agree  
replace c19_stress = . if c19_stress<0
gen c19_stress_agree = inlist(c19_stress,1) if !missing(c19_stress)
* hh lost income (all most some / none)
replace c19_lostincome = . if c19_lostincome<0
gen c19_income_some = inlist(c19_lostincome,1,2,3) if !missing(c19_lostincome)
gen c19_income_most = inlist(c19_lostincome,1,2) if !missing(c19_lostincome)


*** value label of cash arm 
label define treat_cash 0 "IPT-G" 1 "IPT-G+", modify


 
*** Regressions of cash treatment compared to therapy 
*** Regressions
regs c19_getangry c19_argue c19_talkthrough  c19_fear c19_stress_agree c19_income_some, ///
	treatment(treat) rounds(2 3) 
	
*** Q values 
get_qvals c19_getangry c19_argue c19_talkthrough c19_fear c19_stress_agree c19_income_some, rounds(2 3) tname_in("1.treat") tsuff_out("ther") method("simes")
get_qvals c19_getangry c19_argue c19_talkthrough c19_fear c19_stress_agree c19_income_some, rounds(2 3) tname_in("2.treat") tsuff_out("cash") method("simes")
get_qvals c19_getangry c19_argue c19_talkthrough c19_fear c19_stress_agree c19_income_some, rounds(2 3) tname_in("est_diff") tsuff_out("diff") method("simes")


*** unadjusted table
* 12m:
#delimit ;
esttab 	c19_getangry_2 c19_argue_2 c19_talkthrough_2 c19_fear_2 c19_stress_agree_2 c19_income_some_2 
	using "${tables}/results-covid-revised-2.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat 2.treat) 
	stats(c_mean c_sd N pval_ther qval_ther pval_cash qval_cash est_diff pval_diff qval_diff, fmt(3 3 0 3 3 3 3) 
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values" "H0: IPT-G+=0 p-values" "\hspace{20pt}FDR adj. q-values" "Coeff.: IPT-G+ - IPT-G" "H0: IPT-G=IPT-G+ p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 
* 24m: 
#delimit ;
esttab 	c19_getangry_3 c19_argue_3 c19_talkthrough_3 c19_fear_3 c19_stress_agree_3 c19_income_some_3
	using "${tables}/results-covid-revised-3.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat 2.treat) 
	stats(c_mean c_sd N pval_ther qval_ther pval_cash qval_cash est_diff pval_diff qval_diff, fmt(3 3 0 3 3 3 3) 
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values" "H0: IPT-G+=0 p-values" "\hspace{20pt}FDR adj. q-values" "Coeff.: IPT-G+ - IPT-G" "H0: IPT-G=IPT-G+ p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 


